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Abstract. 

I present a simple model of the time dependence of the contact area between 
solid bodies, assuming either a totally uncorrelated surface topography, or a self 
affine surface roughness. Time dependent deformation due to "creep" processes is 
incorporated using a recently proposed model, and produces the time increase of the 
contact area A(t). For an uncorrelated surface topography, A(t) is numerically found to 
be well fitted by expressions of the form [A{oo) — A{t)] ~ {t + to)~'^, where the exponent 
q depends on the normal load as g ~ F^, with /3 close to 0.5. In particular, when the 
contact area is much lower than the nominal area I obtain A{t) /A{0) ~ 1+C ln(t/to+l), 
i.e., a logarithmic time increase of the contact area, in accordance with experimental 
observations. The logarithmic increase for low loads is also obtained analytically in 
this case. For the more realistic case of self affine surfaces, the results are qualitatively 
similar. 
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1. Introduction 

The contact area between solid bodies is in general only a tiny fraction of the apparent 
(or nominal) macroscopic contact area, and is essentially proportional to the normal 
force between the bodies[ll|2]. This fact is at the core of many important characteristics 
of friction phenomena. In particular, it gives a natural explanation of the independence 
of the friction force with the nominal contact area, and the linear increase of friction 
force with normal load. However, contact area, even under static external condition is 
a time dependent quantity, and this produces time dependent effects in friction. For 
instance, the logarithmic time increase of static friction coefficient with the contact 
time, that is well described in many materials [H H], is mostly related to a corresponding 
logarithmic increase in time of the contact area [5]. Direct evidence of this geometric 
aging effect was provided by Dieterich and Kilgore|6]. The time increase of contact 
area is also related with the phenomenon of velocity weakening during sliding, namely 
the fact that the friction force between sliding surfaces can decrease when the relative 
velocity is increased. The velocity weakening effect is extremely important in relation 
with the dynamics of the friction process. In fact, it can be shown that this effect is at 
the base of the process of earthquake generation during the relative motion of tectonic 
plates[7], and also plays an important role in the stick-slip dynamics of many sliding 
systems[I]. These examples highlight the central role that the time increase of contact 
area has in the description of many frictional phenomena. An understanding of the time 
dependence of contact area is thus a very basic aim of any friction theory. 

It is generally accepted that the physical mechanisms by which contact area 
increases in time are associated with plastic phenomena occurring in the materials. For 
our purposes, plastic processes can roughly be classified into two qualitatively different 
groups [U [9]. Rapid plastic effects occur when the imposed deformation conditions 
produce the overpassing of the yield stress of the material. This behavior is typically 
referred to as plastic flow. For lower applied stresses, in cases in which the yield stress is 
not overpassed, there still exists the possibility of thermally activated reacommodations 
in the system that tend to reduce gradually its free energy. These processes are much 
slower than those of plastic flow, and are strongly temperature dependent. They 
are refereed to as creep phenomena. When two solid bodies are pressed together, 
the real contact occurs only in a tiny part of the surface, and the local stresses are 
correspondingly very high. This may produce local plastic flows that increase the 
contact area until the local stresses decay below the material yield stress. After this 
initial stage, the contact area typically increases slowly in time due to creep processes. 
The time increase of contact area in this stage is seen to be logarithmic in time. 

Creep processes (that in other contexts are referred to as aging effects) produce 
in the contacting bodies a tendency to reach progressively more stable conflgurations. 
The meaning of "more stable" here, is that the atomic rearrangements implied by a 
creep process must produce on average a decrease in the total free energy of the system. 
This stabilizing tendency originates in the fact that energy barriers to jump onto a 
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lower energy state are in general lower than those to go to higher energy states, so on 
average creep produces always an energy reduction in the system. Along these lines, the 
phenomenon of contact area increase has been qualitatively described in some idealized 
geometries |10[ |5], and the experimental logarithmic increase law has been justified for 
these cases. Crucial in this analysis is the assumption of a phenomenological creep law 
in which some generic strain rate s is exponentially related to an appropriate stress 
a, i.e., s ~ exp(o"/5'), where S is called the strain rate sensitivity and is temperature 
dependent [5]. The modeling of the phenomenon in more realistic cases is prevented by 
the complication to describe creep in the materials in a sensible and analytically (or even 
numerically) tractable way. As far as I know there is no statistical model that, based 
on well defined microscopic evolution laws, is able to predict the logarithmic increase of 
contact with time. 

Recently, a way to incorporate creep effects in the dynamics of sliding friction 
has been proposed in [HI [12]. In particular this has been applied to models (the 
Burridge-Knopoff|T3] and 01ami-Feder-Christensen|14j models) used to describe seismic 
phenomena. This kind of models describe the friction phenomena occurring at the flat 
(on average) contact surface between two solid bodies that are sheared against each 
other. The kind of modeling proposed in [Til [12] consists in defining some "plastic" 
degrees of freedom, and arguing that they evolve with a tendency to minimize the total 
energy of the system, which is a generic realization of creep processes as described 
in the previous paragraph. It was shown that the proposed relaxation mechanism is 
able to generate realistic sequences of earthquakes, a goal that had not been obtained 
previously without this kind of modification. Also, realistic frictional properties are well 
reproduced using this relaxation mechanism. In particular, a logarithmic increase of the 
static friction coefficient with contact time, and an approximately logarithmic decrease 
of the average friction force as a function of relative velocity has been obtained. 

In view of the standard interpretation of macroscopic friction features in terms of 
contact area, and since the structural relaxation mechanism in [HI [12] is consistent 
with macroscopic friction properties, the question arises if there is a way to use that 
relaxation mechanism to model the time increase of contact area of solid bodies under 
static contact. Such a modeling, if successful, would give further support to the 
structural relaxation mechanism, and would provide an appropriate framework to study 
the phenomenon of time increase of contact area in greater detail than the rather 
qualitative descriptions available up to now. I attempt this particular goal here. In 
order to do this I explain a variation of the models presented in [TTl IT2] that allows 
to define the contact area, and at the same time incorporates the structural relaxation 
mechanism. I show that the main phenomenology associated to the contact area is 
re-obtained, particularly, its logarithmic time increase. 
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2. The Model 

Typically, the surfaces of solid bodies in contact are at the same time rough and elastic 
(or elasto-plastic), and this is important when studying sliding friction. In the context 
of static contact however, some simplifying assumption can be made. 1 will consider the 
case of an elastic surface that is perfectly fiat down to atomic scale in the absence of 
external forces, and an opposing surface that is atomically rough, but strictly rigid. The 
two surfaces are oriented horizontally. The underlying rigid rough surface is described 
by a random variable defined over the plane ^(r). In the numerical simulations the 
values of r will be restricted to lie on a two dimensional square mesh. Two cases will 
be considered separately: one in which the ^(r) are drawn from a unique Gaussian 
distribution independently for each value of r, and a second case in which the ^(r) is 
spatially correlated in order to model a self-afiine surface (see below). 

A realistic modeling of the upper elastic surface should consist in principle in 
determining the equilibrium values of a three dimensional vector displacement field u 
depending upon the two horizontal coordinates, under the action of the surface forces, 
taking into account the elastic response of the surface, measured by an appropriate 
response function. This is what a full contact mechanics calculation aims to. In the 
present case, and already foreseeing the further inclusion of creep effect, some drastic 
simplifications have to be made. 

The simplified description of the elastic surface (sketched in Fig. 1) consists of a 
collection of scalar coordinates u{r), representing the vertical positions of the elastic 
surface at point r, which are coupled via elastic springs connecting nearest neighbor 
mesh points, and elastic interactions with the bulk of the material, defined by a 
set of coordinates M°(r). Note that the only degrees of freedom I allow are vertical 
displacements, and along this direction all the springs act. This modeling of the elastic 
surface is highly simplistic, but at this point this is necessary in order to have a solvable 
model. Note in particular that the consideration of only vertical displacements means 
that we are dealing with a material with "zero Poisson ratio" . Also, it can be seen that 
a localized force onto this surface generates a distortion that decays exponentially with 
distance, whereas the true response of a semi-infinite elastic body is known to decay as 
~ r^^. The limitations and some unrealistic features of the elastic model I am using are 
described in more detail in the last section of the paper. 

The restriction imposed by the contact geometry is that u{r) > ^(r). Given a set 
of values u^{r), the elastic energy of the system is 



where (r, r') stands for pairs of neighbor sites on the numerical mesh. In all cases 
periodic boundary conditions will be used. The equilibrium values of m's are thus found 
by solving the set of equations obtained by minimizing Equation ([1]), namely 




(1) 




(2) 
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(where is the discrete Laplacian operator on the square lattice, lattice parameter 
is taken as unit of length) with the constraint imposed by the contact condition 
'w(r) > The number of points for which u{r) = ^(r) is the number of contact 

points, and I will take this number as a measure of the contact area A in the model. 
The global vertical position of the elastic surface can be measured by the mean value 
of the M°'s. Note also that can be interpreted as a measure of the nominal distance 
between the two bodies and that the normal force Fjv is the sum of the forces of all 
vertical springs, i.e, F/v = ki J^ri'^i''^) ~ 

Conceptually in the same manner as in previous workpTl [12], plastic relaxation 
is incorporated through a time dependence of the values of M°(r). The form of the 
evolution will be obtained from the prescription that the energy of the system E tends 
to be reduced during relaxation. Concretely, I use 



1. e, a standard relaxation equation that tends to reduce the value of E over time as 
much as possible. The time scale for relaxation (controlled by R~^) will be assumed to 
be much larger than the elastic time scale in which the elastic variables u accommodate 
to satisfy Eq. ([2]). Note that the meaning of Eq. ([3]) is that relaxation makes the force 
exerted by the vertical springs tend to an uniform value when t — )■ oo. 

The variable that is assumed we can control is u^. In this respect, note that this 
variable is not modified by the time evolution [Eq. ([3])]. Actually, this is one of the 
reasons to use a conserving dynamics, in which the Laplacian operator is introduced 
in Eq. (|3]), instead of a non- conserving one, in which the Laplacian is absent [15] (in 
connection with this choice, see also the final Section of the paper). To describe a 
possible experimental situation, I first assume (Fig. [11(a)) that is very large, in 
such a way that there is no contact between the two surfaces at any point. Allowing 
infinite time to relax under this condition, the system reaches the uniform state in which 
u°(r) = u'^ and u{r) = everywhere. This is in fact the most relaxed configuration 
since the elastic energy of every spring is zero. I then place at some value in which 
contact occurs at some positions (Fig. [H^b)), and solve Equations ([2]) and ([3]) as a 
function of time. Numerically, the procedure consists in advancing the solution of Eq. 
([3]) for vP by one time step, solve Eq. ([2]) using a standard relaxation algorithm[16] for 
the new values of m, and iterate the process. 

3. Results for uncorrelated roughness 

I first show results in the case in which the variables ^ (r) that describe the roughness of 
the underlying surface are taken independently at each site, from a Gaussian distribution 
of zero mean and unitary variance. Although this case is not very realistic, we will see 
that in addition to the possibility of accurate numerical simulation, it allows for very 
insightful analytical treatment. Results corresponding to i? = are presented in Fig. 

2, where the contact area A (i.e., the number of points in contact) is plotted against 




(3) 
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Figure 1. One dimensional sketch of the modeL (a) A flat and elastic surface 
(defined by the black dots) on top of a rigid rough one (defined by the top of the 
vertical segments), before contact. (b) The situation after contact. Dots [with vertical 
coordinates u(r)] are in equilibrium under the action of vertical and horizontal springs, 
and eventually the force exerted by the rigid rough surface defined by ^(r). The 
coordinates u°(r) evolve in time according to the relaxation equation ([3|). The external 
control parameter is the mean value u^. 



-F/v- Curves for different values of the ratio ki/ko are shown. Two main regimes are 
observed. For low normal load the contact area is essentially proportional to the load, 
whereas if load is too high, we reach a regime of full contact. The crossover between 
partial contact and full contact occurs at a value of F/v that depends on the elastic 
constants of the model. If one of the spring constant dominates over the other the 
crossover value of is proportional to the dominating spring constant. For instance, if 
ko is negligible compared to ki, the elastic surface becomes a collection of independent 
springs (in other contexts this kind of description of an elastic surface is described as 
the Winkler modeljTT]), and the crossover to full contact occurs for a normal load per 
spring of order kia, where a is the typical roughness of the surface, that is taken as 1 
in the present simulations. 

Given the equilibrium configuration for some value of F^, we can set a finite value 
of R and follow the evolution of the contact area in time. In this process it has to 
be taken into account that if we keep fixed, the value of Fn will change in time. 
Since experiments are usually done at constant normal force instead of constant relative 
distance, I implemented a feedback loop in the simulation that allows to keep the value 
of F/v as constant by changing m^. 
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Figure 2. Contact area as a function of normal load, for different values of ki/ko 
without relaxation (full symbols), and the fully relaxed state (open symbols), in which 
all the forces on the ki springs are equal, in a system of 200 x 200 elements. Note 
that the relaxed state has always a contact area larger than unrelaxed states with the 
same normal force. For reference, a dotted line with a slope of 1 is also plotted. 



The asymptotic {t — t- oo) value of the contact area is determined (according to 
Eq. [3]) by the condition that forces on all springs ki are equal [H]. In this situation the 
actual value of ki is irrelevant, and the fully relaxed contact area becomes a function 
of Fjsf/ko (Fig. 1, open symbols). Note that this asymptotic value of the contact area 
is not in general equal to the nominal area, i.e., the system does not reach full contact 
unless F/v overpasses a crossover value that is proportional to k^. In particular this 
means that if we consider the simplest case fco = we obtain the unrealistic result that 
the contact area tends to the nominal value as t — )■ oo. This is the reason that makes 
mandatory the use of a lateral spring fco in the model. In view of standard experimental 
conditions, below I will concentrate in cases in which F/v is such that the contact area is 
only a small fraction of the nominal one. I found that the dependence of the asymptotic 
contact area with normal force at low loads follows a power law A{t oo) ~ F^, with 
an exponent p ~ 0.8. 

The curves in Fig. [2] correspond either to "instantaneous" or to "fully relaxed" 
values of the contact area. As a function of the contact time, the values of the contact 
area must evolve between these two limits. First of all, note from Fig. [2] that the fully 
relaxed contact area is always larger than the unrelaxed configuration with the same 
normal force. This means that contact area will increase in time when relaxation is 
acting, which is the expected result. The full temporal evolution of the contact area 
obtained by solving Eqs. ([2]) and ([3]) is shown in Fig. 3 for different values of the normal 
force Fjv. I have found that all curves of Fig. 3 are very well described by expressions 
of the form 

A(t)M(0) = a - , , (4) 

i.e., they indicate a saturation towards the asymptotic value with the form of a power 
law. The values of a, g, and to depend on Fat. In particular, the dependence of the 
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Figure 3. Time evolution of tlie contact area, for different values of F/v, and two 
values of the ratio fci/fco, for a system of 200x200 sites. The corresponding fitting 
with the expression given in the text (Eq. |4]) is shown by thin continuous lines. The 
exponent q in the fitting function is plotted as a function of Fn in the inset, where it 
is seen that q goes to zero with Fn as a power law (dotted line in the inset has a slope 
of 0.5). 

exponent q on is seen in an inset in Fig. El It is remarkable that q has a dependence 
on Fn of the form q ~ F^ with f3 being close to 0.5. This means in particular that for 
F/v — )■ 0, we can approximate Eq. @ by 

A{t)/A{0) = 1 + C \n{t/to + 1) for Fn ^ (5) 

with C = {a — l)q. It is numerically found that the value of C is not singular in this limit. 
Since experimentally the values of Fj^ are typically tiny compared to those necessary 
to reach full contact, we can say that Eq. (jS]) shows that in general we must expect a 
logarithmic increase of the contact area in time. 

4. Results for Self-AfRne Surfaces 

The results of the previous section were obtained using uncorrelated asperities. 
Numerically, there is no additional complication in trying more realistic distribution 
of surface roughness (although this will preclude analytic treatments as that of the next 
Section). Real surfaces are in fact better described as self affine fractals[T9], and are 
characterized in terms of its Hurst exponent H. This exponent measures the decaying 
in wave vector of the spectral distribution of surface roughness. In this section I present 
numerical results using a self affine surface, and show that the results are qualitatively 
similar to those of the previous section. 

I construct the self affine rough surface using the successive random mean point 
algorithm of Voss |20l |2T] . The surface is characterized by its Hurst exponent H and its 
small scale rms roughness A. The algorithm for the definition of the self- affine rough 
surface proceeds as follows [2Tj: given a mesh of size L x L (for convenience L is chosen 
to be a power of 2) the central point of the mesh is given a value of ^ chosen at random 



Towards a modeling of the time dependence of contact area between solid bodies 



9 




■ 100 

■ 1000 



Figure 4. (a) Topography of one realization of a self affine surface {H = 1/2) in a 
a system of size 256 x 256 sites and ki/ko = 10. Surface height goes form -0.3 to 
0.3 from back to white, (b) Contact surface as a function of load, in the absence of 
relaxation. 

from a Gaussian distribution with zero mean and width i^A, where i = L/\/2 is the 
distance from the center to the corners. This center now becomes a corner of four new 
squares rotated by 45 ° and with a new center-to-corner distance i smaller by a factor 
\/2. The value of ^ at the center of each new square is obtained as the average of ^ 
at the four corners plus a random value chosen from a Gaussian of width i^A. The 
process is iterated down to £ = 1. This algorithm produces a surface that is self-affine 
in the spatial scale between the mesh size and the full system size. A Hurst exponent 
H = 0.5 will be used (tests using other values of H show no qualitative differences). 
I also use the value A = 0.01. An example of the kind of surface obtained by this 
method is presented in Fig. Il](a). Once the underlying rough surface ^(r) is defined in 
this way, the contact with the elastic surface is numerically evaluated exactly by the 
same methods used in the previous Section. The actual contact area in the absence of 
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Figure 5. Contact area for a self affine surface with Hurst exponent H — 1/2, 
using ki/ko — 10. Solid symbols: Contact area as a function of normal load without 
relaxation, for two system sizes. A power law with exponent 0.8 is plotted as a 
reference. Note the independence of the contact area on system size in the region 
away from full contact. Open symbols: corresponding fully relaxed states for the two 
system sizes analyzed. 



relaxation and for different values of F/v can be seen in|4](b). This looks qualitatively 
similar to the results obtained using more realistic elastic surfaces, as that made in 
Ref. [21]. In Fig. |5]I present results for the dependence of the contact area on Fjq in 
the present model. In the absence of relaxation, the dependence can be fitted by an 
expression of the form A ~ F^, with q ~ 0.8, i.e, is a slightly sub-linear dependence. 
Other numerical analysis of this problem, like those in Ref. [21j have obtained a linear 
dependence. I attribute this slight discrepancy to the somewhat artificial description of 
the elastic surface in the present approach, compared with the more precise, truly three 
dimensional description in Ref. [2T] . 

Now I will consider the effect of relaxation. The value of the limit contact area, 
i.e, the contact area at infinite time as a function of the normal force can be seen also 
in Fig. IHl This curve, as in the uncorrelated case, was obtained in a simulation in 
which a constant force on each of the vertical springs is imposed. The results for the 
increase of contact area with time are presented in Fig. |6] (note that the values of -F/v/fco 
correspond, according to Fig. [5] to cases in which the contact area is a small fraction of 
the nominal one). We see that qualitatively the behavior is very similar to the previous 
case. In particular, for the case of very light loads and except for very long times, a 
very good fitting is provided by an expression as that given by Eq. (j5]), with C being 
independent of the external load. 

We can see the actual contact surface at different times in Fig. [71 We observe how 
the increase of contact area involves both an increase of the area of individual contact 
spots, and the appearance of new ones. This trend is very similar to the experimental 
findings of Dieterich and Kilgore|6]. The asymptotic surface of contact (which cannot be 
typically accessed experimentally), shown in the last plot in Fig. [71 displays a uniformly 
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Figure 6. Time evolution of tlie contact area in the presence of relaxation for a self 
affine surface with H ~ 1/2 (a) Un-normalized values for different normal forces and 
system sizes, (b) Results in units of the initial area for two different values of ki/ko. 
The results in (b) are well approximated by a logarithmic time increase with a slope 
that is independent of normal force, and roughly proportional to the value of ki/ko- 
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Figure 7. Time evolution of the contact surface for Fn /ko — 10, ki/ka = 10. System 
size is 256 x 256. The final configuration in the last panel {t — > oo) is obtained from 
an independent simulation, as explained in the text. 
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scattered distribution of contact points. In particular, points that were at contact in 
the first stages of the process may become detached at very long times, due to the stress 
redistribution that occurs due to relaxation. 



5. Analytical Results 

In the limit of very light loads, the contact between surface and substrate occurs only in 
very few points. If we are considering the case of uncorrelated surface roughness, these 
contact point will typically be well separated spatially. This allows for an analytical 
solution of the model in this limit. We will confirm in this way the logarithmic increase 
of contact area in time. The analytical treatment I present turns out to be formally 
similar to that made in the viscoelastic Greenwood- Williamson model [5l [221 [231 Elj . 
fact, the creep phenomena I am modeling through the relaxation mechanism is a kind 
of viscoelastic relaxation [9l [251 EE]. I will come back to this point by the end of the 
paper. 

The strategy to find the solution in this limiting case is to calculate the response 
of the system in an "indentation hardness test" [26] , and then exploit the linearity of 
Eqs. ([2l[3]) to find the full solution. In fact, due to their linearity, Eqs. ([2l[3]) can be 
solved by Fourier decomposition. In the case there is a single contact point between the 
two surfaces (supposed to be the coordinate origin) on which the force is zero for t < 0, 
and takes some constant value /o for t > 0, a direct calculation of the Fourier modes Uq 
gives: 



Uq = fo 



(6) 



where 



_ Rq^kpki 

and the initial condition u*^(r) = at t = (corresponding to a totally relaxed 
configuration before contact) has been used. 

By Fourier inverting this expression we can obtain the time evolution of the 
surface under this indentation condition. Of particular importance will be the time 
evolution of the variable u at the contact point, i.e, u{r = 0,t), and its velocity 
Vo(t) = du{r = 0,t)/dt. The structure of Eq. ([6|) allows to write the solution in 
the form 

u{r = 0,t) = ^G('^,koRt] (8) 







and 




koRt (9) 



where G{xi,X2) is a dimensionless function of two variables, and G'{xi,X2) = dG/dx2- 
In Fig. [H] I show the values of G as a function of X2, for a few different values of Xi. 
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Figure 8. The dimensionless G{xi,X2) function (Eq. ([5])) as a function of X2, for 
different values of xi. Note the logarithmic increase for large X2, independent of the 
value of Xi. Dotted line is the function ooln(a;2) + ai, with ao = 0.04, ai — 0.3. 



We see that for X2 sufficiently large (this means according to Eq. ([H]), for t sufficiently 
large) the function has a logarithmic increase, with a slope that is independent of other 
parameters of the model. The form of the function G is all we need to get a full solution 
of the contact problem on the assumption of few well separated contacts. 

In the contact geometry, the value of u(r = 0,t) must be constant and equal to 
^(r = 0). In order to maintain the contact force equal to /o, we have to adapt m° 
accordingly, namely we must have 

IT = (1°' 

To solve the problem of many contact points (in the assumption that the distance among 
them is large compared to the size of the distortion that any contact exerts onto the 
surface) we ffist need to generalize the form of vo{t) for a constant force /o, to a new 
function v{t) for an arbitrary time dependent form of the contact force f{t). This can 
be done from the previous solution because of the linearity of the equations. The result 
is 

v{t) = drR^G' (^^, koRit - r) j (11) 

To generalize Eq. f llOj) to the case of many contacts, under the action of a total normal 
force -F/v, I use the fact that this normal force has to be distributed among all contacts 
at any time. We obtain 

IT ~ N{t) ^ ^ 

where N{t) is the number of points at contact at time t, m labels the contact points, 
and Vm{t) is expression ( fTTI) calculated with the particular (still unknown) form of fm{t) 
for the contact force at point m. Using Eq. ( ITTi) and the fact that all forces must sum 
up to F/v, we obtain the simple result 

du^ vo{t) 



dt N{t) 



(13) 
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where fo(t) is given in Eq. Namely, the change of the control variable is equal to a 
prescribed function of time, divided by the actual number of contacts at that time. From 
here, and assuming a generic distribution of asperities given by a probability distribution 
P(^), the following equation for the time evolution of the number of contacts can be 
derived: 

. ^ MWM) (14, 

where A^o is the total number of mesh points in the system. Integrating, one obtains 

JN{t=0) P{u^{t')) Jo 

The integral on the left of this equation is explicit, but it is not analytic if the P 
distribution is taken to be a Gaussian. We can replace the Gaussian form of P by an 
exponential to obtain a closed form, i.e., assuming P{^) = 7exp(— 7^), we obtain 

N{t) = N{t = 0)(l + ki vo{t')dt'^ . (16) 

where N{t = 0) is the number of contacts at t = 0, and is given (for exponentially 
distributed asperities) by N(t = 0) = 7/0/^^1- Using the asymptotic form of Vq (Eqs. 
([8]), ([9]), and Fig. [8]), we obtain for large t 



N{t) aok 



-ln{koRt). (17) 



N{t = 0) ko 

with ao — 0.04. As numerically shown by Greenwood and Williamson [23] the 
consideration of a Gaussian distribution of asperities does not alter this result in the 
load range in which 1 ^ N{t) <C Nq. In this way, we obtain analytically a result 
that is nicely compatible with the numerical results: for small loads, the contact area 
increases logarithmically in time, and once scaled with the value of the area at t = 0, 
the result is independent of the precise value of the applied force. Note also that the 
numerical results (Figs. [3]and[6]) show an increase of the logarithm prefactor when ki/ko 
is increased, compatible with Eq. ( [T7|) . The main hypothesis to derive Eq. (fT7|) have 
been the uncorrelated distribution of asperities, and the fact that elastic distortions of 
the surface at the contact points do not influence other contact points, and this means 
that the contact points have to be sufficiently away from each other. I will analyze this 
expression further in the last section of the paper. 



6. Summary and discussion 

In this paper, a simplified model for the time evolution of the contact area between solid 
bodies in contact has been presented. The goals of the model are twofold. On one side, 
for very light loads and in the case of a totally uncorrelated surface roughness, the model 
can be worked out analytically, and it can be shown that a logarithmic dependence of 
the contact area with time emerges. This means that we are able to go all the way from 
a well defined statistical model to its solution, and obtain logarithmic aging. Secondly, 
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numerical simulations can be done in cases where the applied load is not necessarily 
small, and systematic dependences of the contact area in time others than logarithmic 
have been found in this case. Also, the model allows to study more realistic cases 
in which the roughness of the surface is assumed to be correlated. I now focus on a 
discussion of the process of logarithmic increase in time of contact area, that is the most 
directly comparable with available experimental results. 

The analytical results of the previous Section provide the clearest understanding of 
the origin of such a logarithmic increase within the framework of the present model. 
In fact, this logarithmic increase is originated in the form of the surface response 
to a localized and constant applied load. This can be rephrased in the following 
form: If we push the surface of (our model of) an elastic body with a constant force 
at a single point, and look for the deformation this force produces, the indentation 
increases logarithmically in time due to the relaxation processes considered by the 
model. In this respect, I notice that the logarithmic increase is crucially dependent 
on the dimensionality of the surface. For instance, in the case of a line (i.e., the border 
of a half plane) the same relaxation mechanism would produce a displacement that 
grows with time as t^^^ (this dependence appearing when summing up the q modes of 
Eq. ([6]) in a one- dimensional geometry), which is (at least in principle) discernible from 
a logarithmic increase. It thus remains to be seen if there are experimental realizations 
on this confined geometry configuration that allows to test this prediction. 

I now want to discuss in more detail the analytical expression in Eq. f lT7|l . that 
we saw is well verified in the simulations, and compare it with available experimental 
results. One of the most detailed experimental studies of temporal effects in friction 
measurement have been provided by Baumberger and co-workers in a series of papers 
O 123 I2HI |29]. One important parameter they consider, is the coefficient defined as the 
derivative of the static friction coefficient fig with respect to the logarithm of the hold 
time [HI H], namely B = This hold time is the time during which the two surfaces 

are left in rest contact, before the friction force necessary to start sliding is measured. 
This coefficient (that is independent of the time unit chosen) is typically found to have 
a conserved value for a variety of materials, in a range close to 10~^. 

Following Tabor [2], the friction force Fj^ between two solids can be written as 
Ffr = where as is the so-called shear strength of the interface and A the real area 
of contact. Using also the standard expression Ffr = f^s^N, we can write /i^ = Aas/Fjy, 
where we see that the friction coefficient is directly related to the contact area. In 
particular, if we assume that as/F^ takes some constant value, we can write 



This form gives a direct access to the B coefficient as predicted by the model presented 
in this paper. We see first of all (from Eq. [T7|) that B does not depend on the relaxation 



1 d/^s 1 dA 




fig dint A dint 
and since /is is typically of order one, we can roughly write 
„ 1 dA d(A/A(t = 0)) 



(19) 



A dint dint 
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coefficient R. This seems a bit surprising but is not contradictory: although the contact 
area increase in time is produced by a non-zero value of R, this coefficient enters directly 
in setting the time scale (as it is obvious from Eq. ([3])) but the logarithmic derivative 
of contact area is independent of it. Once this fact is recognized, we may ask to what 
extent the value of B in our model can be considered to be constant, independent of 
other parameters, and if this is the case, if this value is in the experimentally observed 
range B ~ 10~^. In this respect, it seems that the answer is negative, since the coefficient 
is directly related to the ratio ki/ko (Eq. ( |T7I) ). which can in principle be set arbitrarily. 
However, I already stressed the fact that the description of our elastic surface is not very 
accurate. In fact, the elastic properties of the surface of an elastic body, that I consider 
isotropic for simplicity, can be characterized by the values of one elastic constant (the 
Young modulus Y, for instance) and are only weakly dependent of a second parameter, 
namely the Poisson ratio i'. The dependence of B on the dimensionless ratio ki/ko is 
an indication that in a hypothetical more accurate description of the elastic body, B 
cannot depend explicitly on the value of Y. Only a dependence on u can exist. Then 
I expect in fact a rather weak dependence of B on the elastic parameters of the body, 
and then a conserved value of B for different materials. Whether this conserved value is 
compatible with the experimental value ~ 10~^ has to be answered once a more realistic 
description of the elastic surface of the body is done. 

Going a step further we must discuss the effect of temperature. The Brechet-Estrin 
analysis [10] predicts that the B coefficient must be temperature dependent, with values 
that increase with increasing temperature. Berthoud et a/. [27] were able to observe 
systematic variations of this coefficient with temperature in different experimental 
situations, that are roughly compatible with the predictions of Brechet and Estrin. 
Typically, an increase of B (up to a factor of roughly ten) was observed when the glass 
temperature of the material was approached. In the present model it is not obvious 
where a temperature dependence can enter that alter the value of the B coefficient. The 
only obvious temperature dependent parameter is the relaxation coefficient R, but we 
have already seen that B is independent of R. 

A possible way out to this situation may involve to consider the possibility of other 
mechanisms of relaxation, in addition to the one considered in Eq. ([3]). For instance, 
in addition to the present mechanism that responds to the fluctuations of the forces on 
the ki springs, we can add a term that is directly dependent of the force itself. This 
would generalize Eq. ([3]) to an equation of the type 

One of the main qualitative difference caused by the inclusion of the last term is 
that under the action of a constant value of the value of goes to zero at very large 
times (contrary to the flnite value reached in the presence of the flrst term alone). In 
this sense the last term represents a "viscous" relaxation in the system, that may have 
a progressively larger effect as the temperature increases towards the glass temperature 
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Figure 9. Depth of an indentation experiment on the model, under a constant load /o, 
as a function of time, for ki/ko = 3, in the presence of creep relaxation and viscoelastic 
relaxation, measured respectively by the coefficients R and R. (see Eq. (PH)) ). Note 
the increase in the slope of the asymptotic logarithmic behavior, as R/ R is increased. 

of the system. This suggests that R/R may be considered to be an increasing function 
of temperature. 

The effect of this term can be seen in Fig. [91 There I show the time dependence 
of the indentation depth for a point contact on the surface of the elastic body for 
increasing values of R/R. These results are obtained by solving analytically the system 
equations for Ug as I did in Section V, using Eq. (!20!) instead of ([3]) , and Fourier inverting 
numerically the result. We see that the logarithmic increase of u in time is conserved 
in the presence of the R term, and the slope increases as R increases. This slope 
measures directly the value of B in the model. It can be easily shown that B doubles 
its value as R/R goes from zero to large values. This means that the consideration 
of alternative relaxation mechanisms can justify a variation of the B coefficient with 
temperature. Whether the change of the B coefficient observed by Berthoud et a/. [27] is 
related to a change of relaxation mechanism or not remains to be investigated further, 
both theoretically and experimentally. 
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